Differential Equations

An ordinary differential equation or ODE is a mathematical equation containing a function or functions of one independent variable and its derivatives. The term ordinary is used in contrast with the term partial differential equation or PDE which involves functions and their partial derivatives with respect to more than one independent variable.

A Mathematical Model is a description of a system using mathematical concepts and language. ODEs and the PDEs are excellent tools for this purpose. Indeed, in real world applications, the functions usually represent physical quantities, the derivatives represent their rates of change, and the equation defines a relationship between the two. Because such relations are extremely common, differential equations play a prominent role in many disciplines including Dynamical Systems, Engineering, Physics, Biology…

In pure mathematics, differential equations are studied from several different perspectives, mostly concerned with their solutions — the set of functions that satisfy the equation. Only the simplest differential equations are solvable by explicit formulas. However, some properties of the solutions of a given differential equation may be determined without finding their exact form.

If an analytical solution is not available, the solution may be numerically approximated. The theory of Dynamical Systems puts emphasis on qualitative analysis of Systems described by differential equations, while many numerical methods have been developed to determine solutions with a given degree of accuracy. In general, partial differential equations are much more difficult to solve analytically than are ordinary differential equations. However, sometimes, separation of variables allows us to transform a PDE in a system of ODEs.

In this work only ODEs were taken into account.

Ordinary Differential Equations (ODEs)

Let F be a given function of x, y, and derivatives of y. Then an equation of the form:

$$F\left (x,y,y',\cdots y^{(n-1)} \right )=y^{(n)}$$

where $y$ is a function of $x$, $y'= \frac{dy}{dx}$ is the first derivative with respect to $x$, and $y^{(n)}=\frac{d^{n} y}{dx^{n}}$ is the nth derivative with respect to $x$, is called an explicit ordinary differential equation of order $n$.

More generally, an implicit ordinary differential equation of order $n$ takes the form:

$$F\left(x, y, y', y'',\ \cdots,\ y^{(n)}\right) = 0$$

An ODE of order $n$ is said to be linear if it is of the form:

$$y^{(n)}(x)+a_{n-1}y^{(n-1)}(x)+\cdots+a_2y''(x)+a_1y'(x)+a_0y(x)=Q(x)$$$$(1)$$

where $a_0$, $a_1$, $...$, $a_{n-1}$ are constants and $Q(x)$ is a function of the independent variable $x$. If $Q(x)=0$, the linear ODE is said to be homogeneous.

In general, an nth-order linear ODE has $n$ linearly independent solutions. Furthermore, any linear combination of linearly independent functions solutions is also a solution, the general solution. The general solution of a non-homogeneous differential equation is obtained by adding the general solution of the associated homogeneous equation with a particular solution of the given equation. The coefficients of the linear combination of the solutions are obtained from the given initial conditions of the problem: $y(0)$, $y’(0)$,$\cdots$, $y^(n-1)(0)$.

Simple theories exist for first-order and second-order ordinary differential equations, and arbitrary ODEs with linear constant coefficients can be solved when they are of certain factorable forms. Methods such as:

• Method of undetermined coefficients
• Integrating factor
• Method of variation of parameters
• Separable differential equation

are usually used. Integral transforms such as the Laplace transform can also be used to solve classes of linear ODEs. This last method was widely discussed during our Dynamical Systems and Control lectures in order to study first order and second order systems responses to some specific inputs, such as step or sinusoidal inputs.

By contrast, ODEs that lack additive solutions are nonlinear, and solving them is far more difficult, as one can rarely represent them by functions in closed form. Instead, exact and analytic solutions of ODEs are in series or integral form that can be solved using methods such as:

• Successive Approximations
• Multiple scale Analysis
• Power series solution of differential equations
• Generalized Fourier series

While there are many general techniques for analytically solving classes of ODEs, the only practical technique for approximating solutions of complicated equations is to use numerical methods. Graphical and numerical methods may approximate solutions of ODEs and yield information that often suffices in the absence of exact analytic solutions. Such methods include:

• Euler method — the most basic method for solving an ODE
• Explicit and implicit methods — implicit methods need to solve an equation at every step
• Backward Euler method — implicit variant of the Euler method
• Trapezoidal rule — second-order implicit method
• Runge–Kutta methods — one of the two main classes of methods for initial-value problems

The most popular of these are the Runge-Kutta methods. A 4th order method (5th order truncation method) was implemented in this work.

In order to use numerical methods to solve a nth-order differential equation, the first step is to transform the differential equation into a system of $n$ differential equations of first order: Let be $Z$ a vector having as components the function $y$ and its first $n-1$ derivatives with respect to $x$:

$$ Z=\begin{bmatrix}y \\ y' \\ y'' \\ \cdots \\ y^{(n-1)}\end{bmatrix} = \begin{bmatrix}z_1 \\ z_2 \\ z_3 \\ \cdots \\ z_n\end{bmatrix}$$

The differential equation (1) is, then transformed into:

$$\begin{cases}z_1'=z_2 \\ z_2'=z_3 \\ \vdots \\ z_{n-1}'=z_n \\ z_n'= Q(x)-a_{n-1}z_n- \cdots - a_1z_2 - a_0z_1 \end{cases}$$

In Dynamical Systems, the independent variable is always the time so, from this point on, we are going to change the notation $x$ to $t$.

Runge-Kutta method

The Runge–Kutta methods are a family of implicit and explicit iterative methods, which are used in temporal discretization for the approximation of solutions of ordinary differential equations. These techniques were developed around 1900 by the German mathematicians C. Runge and M. W. Kutta.

$$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}$$$$k_{1}=f(t_{n},y_{n})$$$$k_{2}=f(t_{n}+c_{2}h,y_{n}+h(a_{21}k_{1}))$$$$k_{3}=f(t_{n}+c_{3}h,y_{n}+h(a_{31}k_{1}+a_{32}k_{2}))$$$$\vdots$$$$k_{s}=f(t_{n}+c_{s}h,y_{n}+h(a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1}))$$

To specify a particular method, one needs to provide the integer s (the number of stages), and the coefficients $a_{ij}$, $b_i$, and $c_i$. The matrix $a_{ij}$ is called the Runge–Kutta matrix, while the $b_i$ and $c_i$ are known as the weights and the nodes.

From the family of the runge-kutta methods the most commonly used are the order 4 and 5 methods, which can be implemented as follows, by inlining the values for $a_{ij}$, $b_i$, and $c_i$ directly in the code:

Second Order Linear System with 1 degree of freedom

In order to study the potential of the numerical integration of an ODE, the second order linear Dynamical System with different responses to different inputs discussed during DSC lectures was used. This Dynamical System can be mathematically modeled by the following ODE:

$$y''(t)+2 \xi w_n y'(t)+w_n^2 y(t)=F(t)$$

where $w_n$ represents the undamped natural frequency, $\xi$ represents the damping ratio and $F(t)$ a forced exterior action.

The solution $y(t)$ of this kind of differential equations is obtained as a sum of the general solution of the homogeneous differential equation, $y_h(t)$, with a particular solution of the complete differential equation, $y_p(t)$:

$$y(t)=y_h(t)+y_p(t)$$

This is a problem of initial conditions, $y(0)$ and $y'(0)$, which allow the determination of the integration constants.

Natural response of the system

The natural response of the system is obtained in the absence of forced exterior actions, in other words, it is the solution of the homogeneous differential equation:

$$y''(t)+2 \xi w_n y'(t)+w_n^2 y(t)=0$$

However, in order to get a response of the Dynamical System, it is necessary to change the initial conditions to the introduction of an initial perturbation to the System which can be modeled by a Dirac impulse. This is a convenient form to apply Laplace transform method to solve analytically the homogenous differential equation:

$$y''(t)+2 \xi w_n y'(t)+w_n^2 y(t)=\delta(t)$$

In these conditions the transfer function is:

$$G(s)=\frac{1}{s^2+2 \xi w_n s+w_n^2}$$

As referred before, the numerical integration of ODEs requires that each equation of degree $n$ is transformed into a system of ODEs of degree 1. This ODE is of degree 2. Thus, we will transform it into a system of 2 ODEs of degree 1.

This ODE is of degree 2. Thus, we will transform it into a system of 2 ODEs of degree 1, as follows:

$$ z = \begin{bmatrix}z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} y(t) \\ y'(t) \end{bmatrix} $$$$ z' = \begin{bmatrix}z_1' \\ z_2' \end{bmatrix} = \begin{bmatrix} z_2 \\ -2\xi w_n z_2 - w_n^{2} z_1 \end{bmatrix}$$$$(2)$$

The function natural is a builder function that takes as arguments $w_n$ and $\xi$ and returns a lambda function representing the system $z'$:

Now let's consider the following initial conditions: $$y_0 = \begin{bmatrix} y(0)\\ y'(0) \end{bmatrix} = \begin{bmatrix} \sqrt{2}\\ \sqrt{2} \end{bmatrix}$$

We will also consider $w_n=2\pi$.

We will focus on studying the properties of $z_2'$, and plotting only it.

Analytic solution

The characteristic equation for the ODE we are studying can be given by:

$$ s^2 + 2\xi w_n s + w_n^2 = 0$$

Applying the solving formula for polynomials, the roots for this equation are given by:

$$ s_{1,2} = -\xi w_n \pm w_n \sqrt{\xi^2 - 1}$$

So, the general solution of this homogenous equation is:

$$ y(t) = C_1 e^{s_1 t} + C_2 e^{s_2 t} $$

Where $C_1$ and $C_2$ are constants to be determined by the initial conditions.

Because this System has different behaviours depending on the values of $\xi$, we will study the System response with regard to different values of $\xi$, namely for:

• $\xi=0$, an undamped System
• $\xi \in{]0,1[} $, an under damped System
• $\xi=1$, a critically damped System
• $\xi>1$, an over damped System

Note that the characteristic equation equals to the denominator of the transfer function $G(s)$ and thus, the poles of the transfer function are the roots of the characteristic equation.

Undamped system

A System is called undamped when $\xi=0$.

We will study the behaviour of $z$ when $\xi=0$ analytically and numerically.

Analytical solution

Since $\xi = 0$ the ODE can be written in the form:

$$ y(t) = C_1 e^{j w_n t} + C_2 e^{-j w_n t} $$$$ y(t) = A_1 \cos{w_n t} + A_2 \sin{w_n t} $$

And its' first derivative:

$$ y'(t) = -w_n A_1 \sin{w_n t} + w_n A_2 \cos{w_n t} $$

Now considering initial conditions $y(0)$:

$$ y(t=0) = y_0 = A_1$$$$ y'(t=0) = y_0' = w_n A_2 $$$$ A_1 = y_0 $$$$ A_2 = \frac{y_0'}{w_n} $$

Replacing $A_1$ and $A_2$ the expression can be given by:

$$ y(t) = y_0 \cos{w_n t} + \frac{y_0'}{w_n} \sin{w_n t} $$

This is implemented as follows:

Numerical solution

By calling the Runge Kutta function an approximation of the ODE can be calculated:

We will plot the difference between the analytic solution and the numeric solution:

We can see that the errors in the numerical approximation are adding up along time.

Notes

• The natural response of the undamped System is a simple harmonic motion with an amplitude of $\sqrt{y(0)^2 + \frac{y'(0)}{w_n}^2}$ . For the numerical example we used $A=\sqrt{2.05}$;

• As stated before, the analytic method may not be possible to use when in case there is no closed form solution, thus the approximation using the numeric method is generally used as it works for any ODE, but it accumulates errors, as shown before;

• In case the transfer function has two complex conjugated poles located in the $Im$ axis of the complex plan (no real part). Then $y_h(t)=L^{-1} (G(s))$ is a linear combination of two complex exponential functions having only pure imaginary symmetric numbers and, applying the Euler’s formula, the expected solution was an harmonic response. The analytical solution and the numerical integration confirmed this prevision.

Underdamped System

A System is underdamped when $0 < \xi < 1$. Lets set $\xi = 0.05$

In this case the roots of the characteristic equation are 2 complex conjugated numbers, that is to say, in the Laplace transform method, the transfer function has two complex conjugated poles in the 3th and 4th quadrants of the complex plane. Then $y_h (t)=L^{-1} (G(s))$ is a linear combination of two complex exponential functions having conjugated exponents. Then common real exponential part multiplies an expression with the same form of the previous case. So the expected solution was an harmonic response but with a decreasing amplitude. The analytical solution and the numerical integration will confirm this prevision.

The analytical solution with the constants calculated from the initial conditions is given by:

$$y(t)=Ae^{-ξw_n t} cos⁡(w_d t-ϕ)$$

where:

• $w_d=w_n \sqrt{1-ξ^2}$ is the damped natural frequency

• $A=\sqrt{\frac{y'(0)+ξw_n y(0)}{w_d}+y(0)^2}$ is a constant calculated from the initial conditions.

• $ϕ=\tan^{-1}\frac{y'(0)+ξw_n y(0)}{w_d y(0)}$ is the phase angle

As this analytical solution shows, the underdamped System has the interesting property of being enveloped by $env(t)=Ae^{-ξw_n t}$, the upper envelope, and $-env(t)$, the lower envelope.

To compare the analytical solution with the numerical one, in the next plot the numerical solution is in blue and the envelope of the analytical one is in gray, which is implemented as follows:

The System is bounded between $[-env(t), env(t)]$ and so $y(t) \to 0$ when $t \to \infty$.

Critically damped System

A System is called critically damped when $\xi=1$

The characteristic equation has 2 real identical roots. The transfer function in the Laplace method has 1 real pole with multiplicity 2. So $G(s)$ is a sum of two partial fractions, one having in the denominator $s-s_{1,2}$ and the other $(s-s_{1,2})^2$. The first one will give rise to a real exponential function and the other to the same real exponential function multiplied by $t$. So, the natural System response will not be harmonic but $y_h(t) \to 0$ as $t \to \infty$.

The analytical solution is:

$$y(t)=(y(0)+(y'(0)+w_n y(0))t)e^{-w_n t}$$

which confirms what was said in the previous paragraph.

The numerical integration was implemented as:

The System rapidly converges to 0.

Overdamped System

A System is called overdamped when $\xi > 1$

In this case the transfer function, $G(s)$, has two distinct real poles with multiplicity 1. Again, the System natural response will not be harmonic, since the solution of the differential equation will be a linear combination of two exponential real functions. Observing the form of the poles %s_1% and %s_2%, both are the sum of two real numbers being one of them common.

The analytical solution will be of the form:

$$y(t)=e^{-ξw_n t}(C_1 e^{w_n \sqrt{ξ^2-1} t}+C_2 e^{-w_n \sqrt{ξ^2-1} t} )$$

and the constants are evaluated from the initial conditions:

$$C_1=\frac{y(0) w_n (ξ+\sqrt{ξ^2-1})+y'(0)}{2w_n \sqrt{ξ^2-1}}$$$$C_2=\frac{-y(0) w_n (ξ+\sqrt{ξ^2-1})-y'(0)}{2w_n \sqrt{ξ^2-1}}$$

For $\xi=2$ the numerical solution was:

This System converges to 0 but it takes longer to do so.

System response to a permanent harmonic input

We will now focus on the behaviour of the System in the presence of a disturbance. This can be done analytically by studying:

$$y''(t)+2 \xi w_n y'(t)+w_n^2 y(t)=F(t)$$

This can be done analytically by adding to the general solution of the homogeneous equation a particular solution of the complete equation:

$$y(t) = y_h(t) + y_p(t)$$

The type of disturbance applied can be divided into the following:

• Harmonic - The applied disturbance is given by a sinusoidal function
• Periodic - The applied disturbance is given by a fourier series
• Transient - Any other type of disturbance

While the harmonic and periodic Systems can be solved analytically or approximated using series, transient Systems may only be solved using convolution of the integrals, Laplace transforms or numeric integration.

The sinusoidal input we are going to consider is of the form $F(t)=f cos(wt)$, where $f$ and $w$ are, respectively, the amplitude and the frequency of the permanent exterior harmonic excitation.

We will considered the undamped System ($\xi=0$), and $f=30$ in order to study the System response to different values of the harmonic excitation input frequency, $w$.

For numerical integration purposes, we define a forced builder function with parameters $w_n$, $\xi$ and f as a multiplier for the force argument which takes a lambda function as follows:

This force function is added to the second member of the last equation, in the first order System of ODEs defined in (2).

We can now study the System response to different values of $w$. First was considered $ w = \frac{5\pi}{6} < w_n = 2 \pi$.

During the lectures on Dynamic Systems, the influence of $w$ was studied using bode plots.

We can see that the disturbance introduced in the System does not change the bounds of the natural function.

Now lets study the interference caused by the function $ w = w_n$:

In this case the period of the natural function and the disturbance functions match ($w = w_n$), and thus the waves are infinitely amplified.

Finally lets consider the case of $ w = w_n - \epsilon = \frac{9}{10}w_n$

This effect is called the beat effect.

Conclusions

ODEs can be solved using a variety of methods however, for practical purposes, numerical methods such as the Runge-Kutta method allows the estimation of the solution of any ODE in reasonable time.

In order to use numeric methods any System of $m$ ODEs of degree $n$ must be reduced to a System of $mn$ ODEs of degree $1$.

We have shown how to do this for a second order ODE.

If the ODE is the mathematical model of a Dynamical System, its response can be studied as a natural response, which in analytic terms, can be defined as $y(t) = y_h(t)$ and the response to a disturbance function is given by $y(t) = y_h(t) + y_p(t)$, where $y_p(t)$ represents the System response to the permanent input.

We solved an ODE both numerically and analytically and have studied its behaviour regarding different settings of a parameter $\xi$.

ODEs are a very powerful Mathematical Model in which the functions can represent physical quantities, the derivatives represent their rates of change, and the equation defines a relationship between the two. Because such relations are extremely common, differential equations play a prominent role in many disciplines.

Bibliography

Vibrações de Sistemas Mecânicos - Apontamentos Teorico-Práticos, José Dias Rodrigues, Faculdade de Engenharia da Universidade do Porto, 2014

Lecture Slides from Dynamical Systems and Control, 2015

Wikipedia

• Ordinary differential equation
• Runge–Kutta methods